/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// in_namelist.c -- variables and functions related to the main input files [Rapaport, p. 33 and pp. 495-498]

#include "dataStructs.h"
#include "in_namelist.h"

// Global variables -- input parameters:
int runId; // job identifier [Rapaport, pp. 500-504, 513]
int newRun; // logical [Rapaport, pp. 500-504]
int keepTetragonal, keepSquare, rdfWat, stepWatRdf, edp, stepEdp, sizeHistEdp, lpp, stepLpp, sizeHistLpp, wpp, stepWpp, sizeHistWpp;
real extPressure, deltaT, extTemperature, tauT, tauP, soluteMass;
real rCutSolute, // solute-solute cutoff
  rCutSoluteElse, 
  rCutLipLip, // cutoff distance between lipid CG sites
  rCutWatWat, /* cutoff distance for ssd-ssd water interactions */
  rNebrShell;
real springRigidity;
int diffusion, /* logical: 1-> evaluate diffusion, 0-> don't do it */
  writeLipLatMotion, stepDiffuse, latDiff, epp, stepEpp, sizeHistEpp,limitWatRdf,sizeHistWatRdf;  
int limitEdp, limitLpp, limitEpp, limitWpp, limitDiff, centerInputStruct, randSeed, insertSolute;
int applyThermostat, applyBarostat, flexBox, stepNebr;
int adjustRegion; // input parameter; logical, used to require modification of the region edges loaded from previous run
int resetTime;    // input parameter; logical, used to require resetting of stepCount and timeNow variables
int removeSystemTranslation; /* logical variable controlling the removal of the net motion of the mass center of the whole system: 
				- if set to 1 (or any other positive integer), activates 'RemoveSystemTranslation()' 
				- if set to 0 (or any negative integer), does nothing */ 
int removeMonolayersTranslation; /* logical variable controlling the removal of the net 'lateral' motion (i.e., in the xy plane) 
				    of the mass centers of each of the two lipid monolayers: 
				    - if set to 1 (or any other positive integer), activates 'RemoveMonolayersTranslation()'
				    - if set to 0 (or any negative integer), does nothing */
int loadStructure, loadVelocities; /* logical */


/* Local (to this module) variables (declared as "static"): */
static const real degrees_IN_radians = ( PI / 180 ); // conversion factor
static char fileNameMD[ 80 ] = "brahms.md"; // name of main input file
static char fileNameAN[ 80 ] = "brahms.an"; // name of analysis input file
static char fileNameLipidFF[ 80 ] = "lipid.ff"; // name of lipid forcefield input file
static char fileNameWaterFF[ 80 ] = "water.ff"; // name of water forcefield input file

/* External data references: */
extern VecR region, regionAdjusted, solInertia;
extern VecI initUcell, initHalfCellWat;
extern real scaleRefSpringLength; 
extern real sigChol, sigAmine, sigPhos, sigEst, sigGly, sigWat; 
extern real epsChol, epsAmine, epsPhos, epsGly, epsWat;
extern real hBondAmineGly, hBondWatWat;
extern real hBondWatEst, hBondWatPhos, hBondWatGly;
extern real hBondAmineWat, hBondAminePhos, hBondAmineEst;
extern real epsTail, epsEst, sigTail, angleRigidity, orientRigidity;
extern real refAngle_cisUnsat, refAngle_PhosGlyEst, refAngle_CholPhosGly;
extern real massChol, massAmine, massPhos, massGly, massEst, massTail, massWat;
extern real inertiaWat, inertiaEst, inertiaGly;
extern real cholCharge, phosCharge, glyDipole, estDipole, watDipole, rangeWatRdf; 
extern int adjustRegion, resetTime, doCheckpoint, nTypes, runId, stepCheckpoint, nebrTabFac;
extern int stepLimit, zConstraint, excludeBondedNebrs;
extern int nSites, nWaters, nLipids, nDOPCsDSPCs, nDOPEs, stepAvg, stepPdb, reCenterBilayer;
extern int writeAreaVol, nSolutes;
extern int limitDiffuseAv, nBuffDiffuse, nValDiffuse, stepDiffuse;
extern int limitLatDiffAv, nBuffLatDiff, nValLatDiff, stepLatDiff;

/************************************************************************************************************
 * ConvertInputParametersToBrahmsUnits() -- converts the value of input parameters from physical units to 
 * brahms units. Brahms units are defined in "definitions.h". Many input parameters are already expressed in
 * brahms units, hence they do not need rescaling.  
 ************************************************************************************************************/
void ConvertInputParametersToBrahmsUnits( FILE *fp ) 
{
  fprintf(fp, "ConvertInputParametersToBrahmsUnits... "); fflush(fp);
  deltaT *= fs_IN_BRAHMS_UNITS;
  tauT *= fs_IN_BRAHMS_UNITS;
  tauP *= fs_IN_BRAHMS_UNITS;
  extTemperature += 273.15; /* convert Celsius (input unit) to Kelvin (Brahms unit) */
  extPressure *= atm_IN_BRAHMS_UNITS; 
  glyDipole *= D_IN_BRAHMS_UNITS;
  estDipole *= D_IN_BRAHMS_UNITS;
  watDipole *= D_IN_BRAHMS_UNITS;
  refAngle_PhosGlyEst  *= degrees_IN_radians; /* convert to radians */
  refAngle_CholPhosGly *= degrees_IN_radians; /* convert to radians */
  refAngle_cisUnsat    *= degrees_IN_radians; /* convert to radians */
  fprintf(fp, "done\n"); fflush(fp);
}

NameList mdList[] = {
  NameR( deltaT ),
  NameI( stepAvg ),
  NameI( stepLimit ),
  NameI( resetTime ),
  NameI( doCheckpoint ),
  NameI( stepCheckpoint ),
  NameI( stepPdb ),
  NameI( runId ),
  NameI( applyThermostat ),
  NameR( extTemperature ),
  NameR( tauT ),
  NameI( applyBarostat ),
  NameR( extPressure ),
  NameR( tauP ),
  NameI( flexBox ),
  NameI( keepTetragonal ),
  NameI( keepSquare ),
  NameI( removeSystemTranslation ),
  NameI( removeMonolayersTranslation ),
  NameR( rCutLipLip ),
  NameR( rCutWatWat ),
  NameR( rCutSolute ),
  NameR( rCutSoluteElse ),
  NameR( rNebrShell ),
  NameI( nebrTabFac ),
  NameI( stepNebr ),
  NameI( nSites),
  NameI( nWaters ),
  NameI( nDOPCsDSPCs ),
  NameI( nDOPEs ),
  NameI( nLipids ),
  NameI( nSolutes ),
  NameI( nTypes ),
  NameR( region ),
  NameI( adjustRegion ),
  NameR( regionAdjusted ),
  NameI( randSeed ),
  NameI( initHalfCellWat ),
  NameI( initUcell ),
  NameI( loadStructure ),
  NameI( loadVelocities ),
  NameI( centerInputStruct), 
  NameI( reCenterBilayer ),
  NameI( zConstraint ),
  NameI( insertSolute ),
};

NameList lipidForcefieldList[] = {
  NameI( excludeBondedNebrs ),
  NameR( sigChol ), 
  NameR( sigAmine ), 
  NameR( sigPhos ), 
  NameR( sigGly ), 
  NameR( sigEst ), 
  NameR( sigTail ), 
  NameR( epsChol ), 
  NameR( epsAmine ), 
  NameR( epsPhos ), 
  NameR( epsGly ), 
  NameR( epsEst ), 
  NameR( epsTail ),
  NameR( hBondWatPhos ),
  NameR( hBondWatGly ),
  NameR( hBondWatEst ),
  NameR( hBondAmineWat ),
  NameR( hBondAminePhos ),
  NameR( hBondAmineGly ),
  NameR( hBondAmineEst ),
  NameR( cholCharge ), 
  NameR( phosCharge ), 
  NameR( glyDipole ),
  NameR( estDipole ), 
  NameR( springRigidity ),
  NameR( scaleRefSpringLength ),
  NameR( angleRigidity ),
  NameR( refAngle_CholPhosGly ),
  NameR( refAngle_PhosGlyEst ),
  NameR( refAngle_cisUnsat ),
  NameR( orientRigidity ),
  NameR( massChol ),
  NameR( massAmine ),
  NameR( massPhos ),
  NameR( massGly ),
  NameR( massEst ),
  NameR( massTail ),
  NameR( inertiaGly ),
  NameR( inertiaEst ),
  NameR( soluteMass ),
  NameR( solInertia ),
};

NameList waterForcefieldList[] = {
  NameR( sigWat ), 
  NameR( epsWat ), 
  NameR( hBondWatWat ), 
  NameR( watDipole ),
  NameR( massWat ),
  NameR( inertiaWat ),
};

NameList analysisList[] = {
  NameI( writeAreaVol ),

  NameI( latDiff ),
  NameI( nBuffLatDiff ),
  NameI( nValLatDiff ),
  NameI( stepLatDiff ),
  NameI( limitLatDiffAv ),
  NameI( writeLipLatMotion ),

  NameI( edp ), 
  NameI( sizeHistEdp ), 
  NameI( stepEdp ),
  NameI( limitEdp ), 

  NameI( epp ), 
  NameI( sizeHistEpp ), 
  NameI( stepEpp ),
  NameI( limitEpp ), 

  NameI( lpp ), 
  NameI( sizeHistLpp ), 
  NameI( stepLpp ),
  NameI( limitLpp ), 

  NameI( wpp ), 
  NameI( sizeHistWpp ), 
  NameI( stepWpp ),
  NameI( limitWpp ), 

  NameI( rdfWat ),
  NameR( rangeWatRdf ),
  NameI( sizeHistWatRdf ),
  NameI( stepWatRdf ),
  NameI( limitWatRdf ),

  NameI( diffusion ),
  NameI( nBuffDiffuse ),
  NameI( nValDiffuse ),
  NameI( stepDiffuse ),
  NameI( limitDiffuseAv ),
};

int ScanNameList( NameList *nameList, char *buffer, int listSize )
{
  int j, k, match, ok;
  char *token;
  FILE *fp;

  if ( ( fp = fopen( buffer, "r" ) ) == 0 ) return ( 0 ); /* return if there's no "brahms.in" file */

  for ( k = 0; k < listSize / sizeof( NameList ); k++ ) {/* scanning input parameters */
    nameList[ k ].vStatus = 0;
  }

  ok = 1;
  while ( 1 ) {
    fgets( buffer, 80, fp ); /* read one line and store it in "buffer" */
    if ( feof( fp ) ) break; 
    token = strtok( buffer, " \t\n" ); /* break the string "buffer" into the substrings delimited by any of " \t\n"  */
    if ( ! token ) break;
    match = 0;
    for ( k = 0; k < listSize / sizeof( NameList ); k++ ) {  /* scanning input parameters */
      if ( strcmp( token, nameList[ k ].vName ) == 0 ) { /* if strings are equal */
        match = 1;
        if ( nameList[ k ].vStatus == 0 ) {
          nameList[ k ].vStatus = 1;
          for ( j = 0; j < nameList[ k ].vLen; j++ ) {
            token = strtok( NULL, ", \t\n" );
            if ( token ) {
              switch ( nameList[ k ].vType ) {
	      case N_I:
		*( ( int * )  ( nameList[ k ].vPtr ) + j ) = atol ( token ); /* convert the string "token" into a number with type "long" */
		break;
	      case N_R:
		*( ( real * ) ( nameList[ k ].vPtr ) + j ) = atof ( token ); /* convert the string "token" into a number with type "double" */
		break;
              }
            } else {
              nameList[ k ].vStatus = 2;
              ok = 0;
            }
          }
          token = strtok( NULL, ", \t\n" );
          if ( token ) {
            nameList[ k ].vStatus = 3;
            ok = 0;
          }
          break;
        } else {
          nameList[ k ].vStatus = 4;
          ok = 0;
        }
      }
    }
    if ( !match ) ok = 0;
  }
  fclose( fp );
  for ( k = 0; k < listSize / sizeof ( NameList ); k++ ) {
    if ( nameList[ k ].vStatus != 1 ) ok = 0;
  }
  return ( ok );
}

/* GetNameLists. Called from "main", reads all the parameters required to
   specify the simulation, topology and analysis from the corresponding input files. */
void GetNameLists( char **argv )
{
  int sizeOfList;
  char buff[ 80 ];

  // read main input file:
  strcpy( buff, fileNameMD ); 
  sizeOfList = sizeof( mdList );
  ScanNameList( mdList, buff, sizeOfList );

  // read water parameters file:
  strcpy( buff, fileNameWaterFF );
  sizeOfList = sizeof( waterForcefieldList );
  ScanNameList( waterForcefieldList, buff, sizeOfList );

  // read lipid forcefield file:
  strcpy( buff, fileNameLipidFF ); 
  sizeOfList = sizeof( lipidForcefieldList );
  ScanNameList( lipidForcefieldList, buff, sizeOfList );

  // read (optional) analysis parameters file:
  strcpy( buff, fileNameAN ); 
  sizeOfList = sizeof( analysisList );
  ScanNameList( analysisList, buff, sizeOfList );
}

void PrintNameList( FILE *fp, NameList *nameList, char *fileName, int listSize )
{
  int j, k;

  fprintf( fp, "Parameters from input file \"%s\":\n", fileName );
  for ( k = 0; k < listSize / sizeof ( NameList ); k++ ) {
    fprintf( fp, "%s\t\t", nameList[ k ].vName );
    if ( strlen( nameList[ k ].vName ) < 8 ) fprintf( fp, "\t" );
    if ( nameList[ k ].vStatus > 0 ) {
      for ( j = 0; j < nameList[ k ].vLen; j++ ) {
        switch (nameList[k].vType) {
	case N_I:
	  fprintf (fp, "%d ", *( ( int * )  ( nameList[ k ].vPtr ) + j ) );
	  break;
	case N_R:
	  fprintf (fp, "%#g ", *( ( real * ) ( nameList[ k ].vPtr ) + j ) );
	  break;
        }
      }
    }
    switch (nameList[k].vStatus) {
    case 0:
      fprintf (fp, "** no data");
      break;
    case 1:
      break;
    case 2:
      fprintf (fp, "** missing data");
      break;
    case 3:
      fprintf (fp, "** extra data");
      break;
    case 4:
      fprintf (fp, "** multiply defined");
      break;
    }
    fprintf (fp, "\n");
  }
  fprintf (fp, "----\n");
  fflush (fp);
}

void PrintNameLists( FILE *fp )
{
  int sizeOfList;

  sizeOfList = sizeof( mdList );
  PrintNameList( fp, mdList, fileNameMD, sizeOfList );

  if ( nLipids ) {
    sizeOfList = sizeof( lipidForcefieldList );
    PrintNameList( fp, lipidForcefieldList, fileNameLipidFF, sizeOfList );
  }

  sizeOfList = sizeof( waterForcefieldList );
  PrintNameList( fp, waterForcefieldList, fileNameWaterFF, sizeOfList );

  if ( fp != stdout ) { 
    sizeOfList = sizeof( analysisList );
    PrintNameList( fp, analysisList, fileNameAN, sizeOfList );
  }
}
